library(RES450)

PRD <- function(x,P){mean(x)*mean(P)/mean(P*x)}
Suits <- function(av,price,nboot=0){
  data <- data.frame(av,price)
  o <- order(data$price)
  data <- data[o,]
  cumav <- 100*cumsum(data$av/1000)/sum(data$av/1000)
  cumprice <- 100*cumsum(data$price/1000)/sum(data$price/1000)
  n = length(cumav)
  dav <- cumav[2:n] + cumav[seq(1,n-1)]
  dav=c(cumav[1],dav)
  dprice <- cumprice[2:n] - cumprice[seq(1,n-1)]
  dprice=c(cumprice[1],dprice)
  L = sum(.5*dav*dprice)
  S =  1 - L/5000
  out = S

  if (nboot>0){
    outse <- array(0,dim=nboot)
    for (i in 1:nboot){
      obs <- sample(1:n,n,replace=T)
      newdata <- data[obs,]
      o <- order(newdata$price)
      newdata <- newdata[o,]
      cumav <- 100*cumsum(newdata$av)/sum(newdata$av)
      cumprice <- 100*cumsum(newdata$price)/sum(newdata$price)
      dav <- cumav[2:n] + cumav[seq(1,n-1)]
      dav=c(cumav[1],dav)
      dprice <- cumprice[2:n] - cumprice[seq(1,n-1)]
      dprice=c(cumprice[1],dprice)
      L = sum(.5*dav*dprice)
      S =  1 - L/5000
      outse[i] = S
    }
    se = sd(outse)
    out <- c(out,se)
  }
  return(out)
}

ExemptSuits <- function(E,m,av,price,g=1){
  taxbase <- pmax(0,m*av - E)
  out <- abs(Suits(taxbase,price))
  return(out)
}

ExemptPRD <- function(E,m,av,price,g=1){
  taxbase <- pmax(0,m*av - E)
  out <- abs(PRD(taxbase/price,price)-g)
  return(out)
}

#### Simulated data
simdata <- function(b,r,n=1000){
  price <- seq(100000,2000000,length=n)
  lnprice <- log(price)
  a = log(r) + (1 - b)*mean(lnprice)
  cat(a,b,"\n")
  lnav <- a + b*lnprice
  av <- exp(lnav)
  ratio <- av/price
  cat("a =",a,"\n")
  print(astats(av,price))
  cat("Suits =",Suits(av,price),"\n")
  for (i in seq(10000,5000000,10000)){
    E <- round(optimize(ExemptSuits,c(i-10000,i),1/median(ratio),av,price)$minimum)
    if (E!=i){break}
  }
  cat("E =",E,"\n")
  out <- list(av,price,E)
  names(out) <- c("av","price","E")
  return(out)
}

csum <- function(x){
  n = length(x)
  cumx <- cumsum(x/n)
  cumx <- cumx/cumx[n]
}
tsum <- function(x){x*sum(tax1)/sum(x)}

sdata <- simdata(.6,n=1000,r=1)
av <- sdata$av
price <- sdata$price
ratio <- av/price
E <- sdata$E
E
lnav <- log(av)
lnprice <- log(price)
PRD(ratio,price)

pdf("./research/assess/exempt/graphs/sim1.pdf")
plot(lnprice,lnav,xlab="Log Sale Price",ylab="Log Assessed Value",lwd=3,
  type="l",ylim=range(lnprice,lnav),main="Log Assessed Value")
lines(lnprice,lnprice,lwd=3,col="red")
legend("topleft",c("Regressive","Statute"),col=c("black","red"),lwd=2)
dev.off()
pdf("./research/assess/exempt/graphs/sim2.pdf")
plot(lnprice,ratio,xlab="Log Sale Price",ylab="Assessment Ratio",lwd=3,
  type="l",main="Assessment Ratio")
abline(h=1,lwd=2,col="red")
legend("topright",c("Regressive","Statute"),col=c("black","red"),lwd=2)
dev.off()


rmed = median(ratio)
### etr when E = 0, 500k
tax1 <- tax1a <- .01*price
tax2a <- tsum(pmax(0, price - 250000))
tax3a <- tsum(pmax(0, price - 500000))
tax1b <- tsum(av/rmed)
tax2b <- tsum(pmax(0,av/rmed - 250000))
tax3b <- tsum(pmax(0,av/rmed - 500000))
taxdata <- data.frame(tax1a,tax2a,tax3a,tax1b,tax2b,tax3b)
etr1a <- tax1a/price
etr2a <- tax2a/price
etr3a <- tax3a/price
etr1b <- tax1b/price
etr2b <- tax2b/price
etr3b <- tax3b/price
etrdata <- data.frame(etr1a,etr2a,etr3a,etr1b,etr2b,etr3b)
head(taxdata)
head(etrdata)

pdf("./research/assess/exempt/graphs/simetr1.pdf")
plot(price/1000,etr1a,type="l",lwd=3,ylim=range(etr1a,etr2a,etr3a),
  xlab="Sale Price ($1000)",ylab="Effective Tax Rate",main="Assessed Value = Sale Price")
lines(price/1000,etr2a,lwd=3,col="red")
lines(price/1000,etr3a,lwd=3,col="blue")
legend("bottomright",c("0","250k","500k"),col=c("black","red","blue"),lwd=2,title="Exemption")
dev.off()

pdf("./research/assess/exempt/graphs/simetr2.pdf")
plot(price/1000,etr1b,type="l",lwd=3,ylim=range(etr1b,etr2b,etr3b),
  xlab="Sale Price ($1000)",ylab="Effective Tax Rate",main="Regressive Assessment")
lines(price/1000,etr2b,lwd=3,col="red")
lines(price/1000,etr3b,lwd=3,col="blue")
legend("bottomright",c("0","250k","500k"),col=c("black","red","blue"),lwd=2,title="Exemption")
dev.off()


rmed = median(ratio)
tax1 <- .01*price
tax2 <- sdata$av
tax3 <- pmax(0,av/rmed - 250000)
tax4 <- pmax(0,av/rmed - 500000)
tax2 <- tsum(tax2)
tax3 <- tsum(tax3)
tax4 <- tsum(tax4)
cprice <- csum(price)
ctax1 <- csum(tax1)
ctax2 <- csum(tax2)
ctax3 <- csum(tax3)
ctax4 <- csum(tax4)

etr1 <- tax1/price
etr2 <- tax2/price
etr3 <- tax3/price
etr4 <- tax4/price



round(Suits(av,price),3)
pdf("./research/assess/exempt/graphs/suits1.pdf")
plot(cprice,ctax1,type="l",lwd=3,xlab="Accumulated Proportion of Total Sale Price",
  ylab="Accumulated Proportion of Tax Payments",ylim=c(0,1),xlim=c(0,1),xaxs="i",yaxs="i")
lines(cprice,ctax2,type="l",lwd=3,ylim=c(0,1),xlim=c(0,1),col="red")
legend("topleft",c("A = P, Exemption = 0:  Area = K","Regressive AV, Exemption = 0:  Area = L"),
  col=c("black","red"),lwd=2)
text(.2,.8,"Suits = 1 - L/K = -0.107")
dev.off()

Suits(tax1,price)
Suits(tax2,price)
Suits(tax3,price)
Suits(tax4,price)

pdf("./research/assess/exempt/graphs/suits2.pdf")
plot(cprice,ctax1,type="l",lwd=3,xlab="Accumulated Proportion of Total Sale Price",
  ylab="Accumulated Proportion of Tax Payments",ylim=c(0,1),xlim=c(0,1),xaxs="i",yaxs="i")
lines(cprice,ctax2,type="l",lwd=3,ylim=c(0,1),xlim=c(0,1),col="red")
lines(cprice,ctax3,type="l",lwd=3,ylim=c(0,1),xlim=c(0,1),col="blue")
lines(cprice,ctax4,type="l",lwd=3,ylim=c(0,1),xlim=c(0,1),col="green")
legend("topleft",c("A = P, E = 0:  Suits = 0","Regressive AV, E = 0:  Suits = -0.107",
  "Regressive AV, E = $250k:  Suits = -0.042","Regressive AV, E = $500k:  Suits = 0.060"),
  col=c("black","red","blue","green"),lwd=2)
dev.off()



rmed = median(ratio)
tax1 <- .01*price
tax2 <- tsum(sdata$av)
tax3 <- tsum(pmax(0,av/rmed - sdata$E))
cprice <- csum(price)
ctax1 <- csum(tax1)
ctax2 <- csum(tax2)
ctax3 <- csum(tax3)
etr1 <- tax1/price
etr2 <- tax2/price
etr3 <- tax3/price

pdf("./research/assess/exempt/graphs/suits3.pdf")
plot(cprice,ctax1,type="l",lwd=3,xlab="Accumulated Proportion of Total Sale Price",
  ylab="Accumulated Proportion of Tax Payments",ylim=c(0,1),xlim=c(0,1),xaxs="i",yaxs="i")
lines(cprice,ctax3,type="l",lwd=3,ylim=c(0,1),xlim=c(0,1),col="red")
legend("topleft",c("A = P, E = 0:  Suits = 0","Regressive AV, E = $363k:  Suits = 0"),
  col=c("black","red","blue"),lwd=2)
dev.off()

pdf("./research/assess/exempt/graphs/suits4.pdf")
plot(price/1000,etr1,type="l",lwd=3,xlab="Sale Price ($1000)",
  ylab="Effective Tax Rate",ylim=range(etr1,etr3))
lines(price/1000,etr3,type="l",lwd=3,ylim=c(0,1),xlim=c(0,1),col="red")
legend("right",c("A = P, E = 0:  Suits = 0","Regressive AV, E = $363k:  Suits = 0"),
  col=c("black","red","blue"),lwd=2)
dev.off()

